5. 行列の基本変形と不変量

Pythonで学ぶ線形代数学

基本行列による単位立方体の変化

code: elementary_vp.py

from vpython import *

import numpy as np

o = vec(0, 0, 0)

x, y, z = vec(1, 0, 0), vec(0, 1, 0), vec(0, 0, 1)

X, Y, Z = o, y, z, y+z, o, z, x, z+x, o, x, y, x+y

box(pos=(x+y+z)/2)

for v in x, y, z:

curve(pos=-v, 3*v, color=v)

def T(A, u):

return vec(*np.dot(A, (u.x, u.y, u.z)))

E = 1, 0, 0], 0, 1, 0, [0, 0, 1

E1 = 1, 2, 0], 0, 1, 0, [0, 0, 1

E2 = 0, 1, 0], 1, 0, 0, [0, 0, 1

E3 = 2, 0, 0], 0, 1, 0, [0, 0, 1

for u, U in (x, X), (y, Y), (z, Z):

for v in U:

A = E

curve(pos=T(A, v), T(A, u+v), color=u)

ソースコード: elementary_vp.py

https://gyazo.com/0aabb55546f713f1ac2908bbe60d7b34

https://gyazo.com/62cd38963f407b80e6cf8c8b9ee89037

https://gyazo.com/f4cd27b6800a25a22a12ad30db030830

基本行列

code: elementary_sp.py

from sympy import Matrix, var

var('x y a11 a12 a13 a21 a22 a23 a31 a32 a33')

E1 = Matrix(1, x, 0], 0, 1, 0, [0, 0, 1)

E2 = Matrix(0, 1, 0], 1, 0, 0, [0, 0, 1)

E3 = Matrix(1, 0, 0], 0, y, 0, [0, 0, 1)

A = Matrix(a11, a12, a13], a21, a22, a23, [a31, a32, a33)

ソースコード: elementary_sp.py

対話モードでの実行例（基本行列を左から掛けた場合））

code: pyhon

>> E1 * A

Matrix([

a11 + a21*x, a12 + a22*x, a13 + a23*x,

a21, a22, a23,

a31, a32, a33])

>> E2 * A

Matrix([

a21, a22, a23,

a11, a12, a13,

a31, a32, a33])

>> E3 * A

Matrix([

a11, a12, a13,

a21*y, a22*y, a23*y,

a31, a32, a33])

対話モードでの実行例（基本行列を右から掛けた場合））

code: python

>> A * E1

Matrix([

a11, a11*x + a12, a13,

a21, a21*x + a22, a23,

a31, a31*x + a32, a33])

>> A * E2

Matrix([

a12, a11, a13,

a22, a21, a23,

a32, a31, a33])

>> A * E3

Matrix([

a11, a12*y, a13,

a21, a22*y, a23,

a31, a32*y, a33])

階数の問題生成プログラム

code: prob_rank.py

from numpy.random import seed, choice, permutation

from sympy import Matrix

def f(P, m1, m2, n):

if n > min(m1, m2):

return Matrix(choice(P, (m1, m2)))

else:

while True:

X, Y = choice(P, (m1, n)), choice(P, (n, m2))

A = Matrix(X.dot(Y))

if A.rank() == n:

return A

m1, m2 = 3, 4

seed(2020)

for i in permutation(max(m1, m2)):

print(f(-3, -2, -1, 1, 2, 3, m1, m2, i+1))

ソースコード: prob_rank.py

実行例

code: python

Matrix(4, 5, -5, 4], 8, 1, -17, 16, [0, -3, -3, 0)

Matrix(-10, 4, -8, -8], 7, -2, 8, 4, [0, -1, -3, 2)

Matrix(-2, 2, -2, 1], -1, -3, -2, -2, [-1, -1, -2, -3)

Matrix(3, -1, -3, -1], 9, -3, -9, -3, [-3, 1, 3, 1)

定義に従って行列式を計算する

code: determinant.py

from functools import reduce

def P(n):

if n == 1:

return [(0, 1)]

else:

Q = []

for p, s in P(n-1):

Q.append((p + n-1, s))

for i in range(n-1):

q = p + n-1

qi, q-1 = q-1, qi

Q.append((q, -1*s))

return Q

def prod(L):

return reduce(lambda x, y: x*y, L)

def det(A):

n = len(A)

a = sum([s*prod([Ai[pi] for i in range(n)]) for p, s in P(n)])

return a

if __name__ == '__main__':

A = 1, 2], [2, 3

B = 1, 2], [2, 4

C = 1, 2, 3], 2, 3, 4, [3, 4, 5

D = 1, 2, 3], 2, 3, 1, [3, 1, 2

print(det(A), det(B), det(C), det(D))

ソースコード: determinant.py

実行結果

code: python

-1 0 0 -18

問5.5の解答例($ n=5)の場合

code: random_matrix.py

from numpy.random import seed, randint

from numpy.linalg import det

def f(n):

S = 0

for i in range(10000):

D = det(randint(0, 10, (n, n)))

if abs(D) < 1.0e-10:

print(D)

f(5)

ソースコード: random_matrix.py

実行例

code: python

0.0

7.087663789206989e-13

0.0

-7.505107646466066e-13

-1.3429257705865887e-12

-1.0924594562311528e-12

-1.3522516439934362e-12

6.03961325396085e-13

非正則行列でも、行列式の値が正確に0にならないことがあることに注意

行列式（逆行列）の問題生成プログラム

code: prob_det.py

from numpy.random import seed, choice, permutation

from sympy import Matrix

def f(P, m, p):

while True:

A = Matrix(choice(P, (m, m)))

if p == 0:

if A.det() == 0:

return A

elif A.det() != 0:

return A

m = 3

seed(2020)

for p in permutation(2):

print(f(-3, -2, -1, 1, 2, 3, m, p))

ソースコード: prob_det.py

実行例

code: python

Matrix(-3, 1, 1], 1, 3, 1, [-3, 3, -3)

Matrix(-2, 1, -1], -3, 2, 3, [3, -2, -3)

連立方程式の問題生成プログラム

code: prob_eqn.py

from numpy.random import seed, choice, shuffle

from sympy import Matrix, latex, solve, zeros

from sympy.abc import x, y, z

def f(P, m, n):

while True:

A = Matrix(choice(P, (3, 4)))

if A:, :3.rank() == m and A.rank() == n:

break

A, b = A:, :3, A:, 3

u = Matrix(x], y, [z)

print(f'{latex(A)}{latex(u)}={latex(b)}')

print(solve(A*u - b, x, y, z))

seed(2020)

m, n = 2, 2

f(range(2, 10), m, n)

ソースコード: prob_eqn.py

実行例

code: python

\left\begin{matrix}7 & 8 & 8\\5 & 8 & 6\\8 & 8 & 9\end{matrix}\right\left\begin{matrix}x\\y\\z\end{matrix}\right=\left\begin{matrix}3\\5\\2\end{matrix}\right

{y: 5/4 - z/8, x: -z - 1}

余因子行列を用いた逆行列の計算

code: inv.py

from numpy import array, linalg, random

n = 5

A = random.randint(0, 10, (n, n))

K = [j for j in range(n) if j != i for i in range(n)]

B = array((-1) ** (i+j) * linalg.det(A[Ki, :][:, K[j)

for i in range(n)] for j in range(n)])

print(A.dot(B/linalg.det(A)))

ソースコード: inv.py

実行結果

code: python

[[ 1.00000000e+00 -2.66453526e-15 -2.66453526e-15 1.59872116e-14

-1.15463195e-14]

[-5.55111512e-16 1.00000000e+00 -1.33226763e-15 4.44089210e-15

-2.22044605e-15]

[-4.44089210e-16 -6.21724894e-15 1.00000000e+00 2.13162821e-14

-1.95399252e-14]

[-1.77635684e-15 -5.32907052e-15 -1.66533454e-15 1.00000000e+00

-1.77635684e-14]

[-1.77635684e-15 -6.21724894e-15 -1.99840144e-15 1.42108547e-14

1.00000000e+00]]

単位行列になっている（非対角要素にわずかな誤差を含む）